#include "dem_reader.h"

DemReader::DemReader() : 
    ncols(0), nrows(0), 
    xllcorner(0.0), yllcorner(0.0), cellsize(0.0), 
    dataLoaded(false) {
}

DemReader::~DemReader() {
    // 清理资源
    for (auto& row : demData) {
        row.clear();
    }
    demData.clear();
}

bool DemReader::loadDemData(const std::string& filePath) {
    std::ifstream file(filePath);
    if (!file.is_open()) {
        std::cerr << "无法打开DEM数据文件: " << filePath << std::endl;
        return false;
    }

    // 清理之前的数据
    demData.clear();
    dataLoaded = false;

    std::string line;
    std::string keyword;

    // 读取头部信息
    while (std::getline(file, line)) {
        std::istringstream iss(line);
        iss >> keyword;
        
        if (keyword == "ncols") {
            iss >> ncols;
        } else if (keyword == "nrows") {
            iss >> nrows;
        } else if (keyword == "xllcorner") {
            iss >> xllcorner;
        } else if (keyword == "yllcorner") {
            iss >> yllcorner;
        } else if (keyword == "cellsize") {
            iss >> cellsize;
        } else if (keyword == "NODATA_value") {
            // 读取NODATA_value，这里我们知道它是"x"
            std::string nodata;
            iss >> nodata;
        } else {
            // 如果不是头部关键字，回退文件指针，准备读取数据
            file.seekg(-static_cast<std::streampos>(line.length() + 1), std::ios::cur);
            break;
        }
    }

    // 读取高程数据
    demData.resize(nrows);
    
    // 由于文件很大，可能需要优化内存使用
    std::cout << "警告: DEM数据非常大 (" << ncols << "x" << nrows << "), 尝试使用更高效的加载方式..." << std::endl;
    
    for (int i = 0; i < nrows; ++i) {
        if (!std::getline(file, line)) {
            std::cerr << "读取DEM数据时出错，行: " << i << ", 已读取" << i << "行" << std::endl;
            return false;
        }
        
        std::istringstream iss(line);
        
        // 预先分配空间
        demData[i].clear();
        demData[i].reserve(ncols);
        
        std::string value;
        while (iss >> value) {
            // 检查是否为无效数据点
            if (value == "x") {
                demData[i].push_back(-9999.0);  // 用-9999表示无效数据
            } else {
                try {
                    demData[i].push_back(std::stod(value));
                } catch (const std::invalid_argument& ia) {
                    std::cerr << "解析DEM数据时出错，行: " << i << ", 值: " << value << ", 设为无效值" << std::endl;
                    demData[i].push_back(-9999.0);
                }
            }
        }
        
        // 检查数据点数量是否符合预期
        // 将ncols转换为size_t类型以避免符号比较警告
        size_t expectedCols = static_cast<size_t>(ncols);
        if (demData[i].size() != expectedCols) {
            std::cerr << "警告: 第" << i << "行的数据点数量 (" << demData[i].size() 
                      << ") 与预期不同 (" << expectedCols << "), 尝试继续..." << std::endl;
            
            // 如果数据点不足，填充-9999
            if (demData[i].size() < expectedCols) {
                demData[i].resize(expectedCols, -9999.0);
            }
        }
        
        // 每100行显示一次进度
        if (i % 100 == 0) {
            std::cout << "已加载" << i << "行数据..." << std::endl;
        }
    }

    file.close();
    dataLoaded = true;
    
    std::cout << "成功加载DEM数据: " << ncols << "x" << nrows << " 网格，左下角坐标(" << xllcorner << ", " << yllcorner << ")，单元格大小: " << cellsize << std::endl;
    
    return true;
}

bool DemReader::latLonToGrid(double longitude, double latitude, int& col, int& row) {
    // 计算列号和行号
    // 列号 = (经度 - 左下角经度) / 单元格大小
    col = static_cast<int>((longitude - xllcorner) / cellsize);
    // 行号 = (左下角纬度 + 行数*单元格大小 - 纬度) / 单元格大小
    // 注意这里行号的计算是倒序的，因为DEM数据通常从上到下存储
    row = static_cast<int>((yllcorner + nrows * cellsize - latitude) / cellsize);
    
    // 检查坐标是否在数据范围内
    // 注意需要-1，因为我们需要至少两个点进行双线性插值
    if (col < 0 || col >= ncols - 1 || row < 0 || row >= nrows - 1) {
        std::cerr << "坐标(" << longitude << ", " << latitude << ") 不在DEM数据范围内" << std::endl;
        return false;
    }
    
    return true;
}

double DemReader::bilinearInterpolation(double x, double y) {
    int col, row;
    if (!latLonToGrid(x, y, col, row)) {
        return 0.0;
    }
    
    // 计算在单元格中的相对位置
    // xRatio和yRatio表示目标点在网格单元内的相对坐标(0-1之间)
    double xRatio = (x - (xllcorner + col * cellsize)) / cellsize;
    double yRatio = (y - (yllcorner + (nrows - row - 1) * cellsize)) / cellsize;
    
    // 获取四个角点的高程值
    // z1: 左下角, z2: 右下角, z3: 左上角, z4: 右上角
    double z1 = demData[row][col];
    double z2 = demData[row][col + 1];
    double z3 = demData[row + 1][col];
    double z4 = demData[row + 1][col + 1];
    
    // 检查是否有无效数据点
    // 无效数据点在文件中标记为'x'，内部存储为-9999.0
    int validPoints = 0;
    double sumZ = 0.0;
    
    if (z1 != -9999.0) { validPoints++; sumZ += z1; }
    if (z2 != -9999.0) { validPoints++; sumZ += z2; }
    if (z3 != -9999.0) { validPoints++; sumZ += z3; }
    if (z4 != -9999.0) { validPoints++; sumZ += z4; }
    
    // 如果所有点都是无效的，返回0
    if (validPoints == 0) {
        std::cerr << "警告: 查询点周围没有有效高程数据" << std::endl;
        return 0.0;
    }
    
    // 如果部分点无效，使用有效点的平均值
    if (validPoints < 4) {
        std::cerr << "警告: 查询点周围部分数据无效，使用平均值" << std::endl;
        return sumZ / validPoints;
    }
    
    // 双线性插值计算
    // 双线性插值的公式：
    // z = (1-x) * (1-y) * z1 + x*(1-y)*z2 + (1-x)*y*z3 + x*y*z4
    // 其中x和y是相对位置，z1-z4是四个角点的高程
    // 这个公式实际上是在两个方向上进行线性插值，先横向，再纵向
    double z = (1 - xRatio) * (1 - yRatio) * z1 + 
               xRatio * (1 - yRatio) * z2 + 
               (1 - xRatio) * yRatio * z3 + 
               xRatio * yRatio * z4;
    
    return z;
}

double DemReader::calculateElevation(double longitude, double latitude) {
    if (!dataLoaded || demData.empty()) {
        std::cerr << "DEM数据未加载" << std::endl;
        return 0.0;
    }
    
    // 使用双线性插值计算高程
    return bilinearInterpolation(longitude, latitude);
}